class: title-slide # ER014 - Data Science & Strategy for Business ## PVA2 ### Teil 3: Stochastische Regressionsanalyse <br> <br> <br> <br> <br> <br> <br> ### FS 2025 <br> ### Prof. Dr. Jörg Schoder .mycontacts[
@FFHS-EconomicResearch
@jfschoder ] --- layout: true <div class="my-footer"></div> <div style="position: absolute;left:400px;bottom:10px;font-size:9px">
Prof. Dr. Jörg Schoder</div> --- name: agenda class: left .blockquote[Agenda] ## Stochastische Regressionsanalyse * Intro * Regressionsanalyse mit Resampling-Methode * Regressionsanalyse im traditionellen Ansatz * Regressionsdiagnostik im traditionellen Ansatz --- class: left .blockquote[Intro] ## Grundidee stochastische Regressionsanalyse .panelset[ .panel[.panel-name[Stichprobe1] <img src="data:image/png;base64,#03_InferenceRegression_files/figure-html/unnamed-chunk-1-1.png" width="55%" style="display: block; margin: auto;" /> ] .panel[.panel-name[Stichprobe 2] <img src="data:image/png;base64,#03_InferenceRegression_files/figure-html/unnamed-chunk-2-1.png" width="55%" style="display: block; margin: auto;" /> ] .panel[.panel-name[Code Sampling] ``` r library(tidyverse) ## Import data ----- my_in_file <- "autos_(StockerUIBK)_20240414.csv" tbl_autos <- read_csv2(xfun::from_root("data","raw",my_in_file)) p <- tbl_autos %>% ggplot(aes(x=Alter,y=Preis)) + geom_point() ## Stichproben ziehen ------ ### Stichprobe 1 ----- set.seed(23) tbl_autos_sub1 <- tbl_autos %>% sample_n(size=7,replace=FALSE) ### "Stichprobe" 2 ------- tbl_autos_sub2 <- tbl_autos %>% slice(c(13,18,21,25)) ``` ] .panel[.panel-name[Code Diagramm] ``` r ## Streudiagramm ---- ### Regressionsgerade gesamt ----- p <- p + geom_smooth(method = "lm", se = FALSE) ### Stichproben ergänzen ----- p + geom_point(data=tbl_autos_sub1, color="#502479") + geom_smooth(data=tbl_autos_sub1, color="#502479",method = "lm", se = FALSE) + geom_point(data=tbl_autos_sub2, color="#d50006") + geom_smooth(data=tbl_autos_sub2, color="#d50006",method = "lm", se = FALSE) ``` ] ] .quelle[Eigene Darstellung.] --- class: left .blockquote[Intro] ## Wiederholung PVA1: mögliche Punktschätzer <table class="table" style="font-size: 16px; margin-left: auto; margin-right: auto;"> <caption style="font-size: initial !important;"> <span id="tab:table-ch8"></span>Mögliche Punktschätzer auf Basis von Stichproben </caption> <thead> <tr> <th style="text-align:right;"> Fall </th> <th style="text-align:left;"> Parameter der Grundgesamtheit </th> <th style="text-align:left;"> Notation </th> <th style="text-align:left;"> Punktschätzung </th> <th style="text-align:left;"> Symbol(e) </th> </tr> </thead> <tbody> <tr> <td style="text-align:right;width: 0.5in; "> 1 </td> <td style="text-align:left;width: 1.2in; "> Anteil in Grundgesamtheit </td> <td style="text-align:left;width: 0.8in; "> <span class="math inline">\(p\)</span> </td> <td style="text-align:left;width: 1.5in; "> Anteil in Stichprobe </td> <td style="text-align:left;width: 0.6in; "> <span class="math inline">\(\widehat{p}\)</span> </td> </tr> <tr> <td style="text-align:right;width: 0.5in; "> 2 </td> <td style="text-align:left;width: 1.2in; "> Mittelwert der Grundgesamtheit </td> <td style="text-align:left;width: 0.8in; "> <span class="math inline">\(\mu\)</span> </td> <td style="text-align:left;width: 1.5in; "> Stichprobenmittelwert </td> <td style="text-align:left;width: 0.6in; "> <span class="math inline">\(\overline{x}\)</span> oder <span class="math inline">\(\widehat{\mu}\)</span> </td> </tr> <tr> <td style="text-align:right;width: 0.5in; "> 3 </td> <td style="text-align:left;width: 1.2in; "> Differenz von Anteilen einer Grundgesamtheit </td> <td style="text-align:left;width: 0.8in; "> <span class="math inline">\(p_1 - p_2\)</span> </td> <td style="text-align:left;width: 1.5in; "> Differenz von Anteilen einer Stichprobe </td> <td style="text-align:left;width: 0.6in; "> <span class="math inline">\(\widehat{p}_1 - \widehat{p}_2\)</span> </td> </tr> <tr> <td style="text-align:right;width: 0.5in; "> 4 </td> <td style="text-align:left;width: 1.2in; "> Differenz von Mittelwerten der Grundgesamtheit </td> <td style="text-align:left;width: 0.8in; "> <span class="math inline">\(\mu_1 - \mu_2\)</span> </td> <td style="text-align:left;width: 1.5in; "> Differenz von Stichprobenmittelwerten </td> <td style="text-align:left;width: 0.6in; "> <span class="math inline">\(\overline{x}_1 - \overline{x}_2\)</span> oder <span class="math inline">\(\widehat{\mu}_1 - \widehat{\mu}_2\)</span> </td> </tr> <tr> <td style="text-align:right;width: 0.5in; "> 5 </td> <td style="text-align:left;width: 1.2in; "> Empirischer Regressionskoeffizient (Grundgesamtheit) </td> <td style="text-align:left;width: 0.8in; "> <span class="math inline">\(\beta_1\)</span> </td> <td style="text-align:left;width: 1.5in; "> Angepasster ("fitted") Regressionskoeffizient (Stichprobe) </td> <td style="text-align:left;width: 0.6in; "> <span class="math inline">\(b_1\)</span> oder <span class="math inline">\(\widehat{\beta}_1\)</span> </td> </tr> </tbody> </table> .quelle[vgl. <a name=cite-ismay_statistical_2024></a>([Ismay and Kim, 2024](https://moderndive.com/index.html))] --- class: inverse, center, middle ## Regressionsanalyse mit der Resampling-Methode .blockquote[Bootstrapping-Verteilung und Konfidenzintervall] .blockquote[Signikanz des Steigungsparameters] --- class: left .blockquote[Bootstrapping-Verteilung und Konfidenzintervall] ## Punktschätzer und Bootstrapping-Verteilung mit der **infer**-Pipe .panelset[ .panel[.panel-name[Koeffizient] * Modell (**ein** numerischer Regressor): `$$\begin{eqnarray} \widehat{\mbox{Preis}}&=&b_0 + b_1\cdot \mbox{Alter}\\ \end{eqnarray}$$` * **infer**-Pipe zur Ermittlung des Steigunskoeffizienten ``` r library(infer) obs_slope <- tbl_autos %>% specify(Preis ~ Alter) %>% calculate(stat = "slope") ``` ``` ## Response: Preis (numeric) ## Explanatory: Alter (numeric) ## # A tibble: 1 × 1 ## stat ## <dbl> ## 1 -2636. ``` ] .panel[.panel-name[Konfidenzintervall] * Bootstrapping-Verteilung des Koeffizienten ``` r ### Bootstrapping-Verteilung des Koeffizienten ---- boot_dist <-tbl_autos %>% specify(Preis~Alter) %>% generate(reps = 1000, type = "bootstrap") %>% calculate(stat = "slope") ``` ] .panel[.panel-name[Diagramm] ``` r visualize(boot_dist) ``` <img src="data:image/png;base64,#03_InferenceRegression_files/figure-html/unnamed-chunk-8-1.png" width="45%" style="display: block; margin: auto;" /> ] ] --- class: left .blockquote[Bootstrapping-Verteilung und Konfidenzintervall] ## Konfidenzintervalle im Vergleich .panelset[ .panel[.panel-name[Perzentilmethode] ``` r percentile_ci <- boot_dist %>% get_confidence_interval(type = "percentile", level = 0.95) percentile_ci ``` ``` ## # A tibble: 1 × 2 ## lower_ci upper_ci ## <dbl> <dbl> ## 1 -2982. -2285. ``` ] .panel[.panel-name[Standardfehler-Methode] ``` r se_ci <- boot_dist %>% get_confidence_interval(type = "se", point_estimate = obs_slope, level = 0.95) se_ci ``` ``` ## # A tibble: 1 × 2 ## lower_ci upper_ci ## <dbl> <dbl> ## 1 -2995. -2276. ``` ] .panel[.panel-name[Diagramm] ``` r visualize(boot_dist) + shade_confidence_interval(endpoints = percentile_ci, fill = NULL, linetype = "solid",color = "dodgerblue2") + shade_confidence_interval(endpoints = se_ci, fill = NULL, linetype = "dashed",color = "dodgerblue3") ``` <img src="data:image/png;base64,#03_InferenceRegression_files/figure-html/unnamed-chunk-11-1.png" width="33%" style="display: block; margin: auto;" /> ] ] ??? sehr ähnlich, keines enthält Null. Legt signifikant negativen Koeffizienten nahe --- class: left .blockquote[Signikanz des Steigungsparameters] ## Konfidenzintervalle im Vergleich .panelset[ .panel[.panel-name[Null-Verteilung] ``` r null_dist <-tbl_autos %>% specify(Preis~Alter) %>% hypothesize(null="independence") %>% #Nullhypothese als zusätzliches Element der Pipe generate(reps = 1000, type = "permute") %>% # permute statt bootstrap calculate(stat = "slope") null_dist %>% get_p_value(obs_stat = obs_slope,direction = "less") ``` ``` ## Warning: Please be cautious in reporting a p-value of 0. This result is an approximation ## based on the number of `reps` chosen in the `generate()` step. ## ℹ See `get_p_value()` (`?infer::get_p_value()`) for more information. ``` ``` ## # A tibble: 1 × 1 ## p_value ## <dbl> ## 1 0 ``` ] .panel[.panel-name[Diagramm] ``` r null_dist %>% visualize() + # Verteilung um Null! - wg. $H_0:~\beta_1=0$ shade_p_value(obs_slope, direction = "less") ``` <img src="data:image/png;base64,#03_InferenceRegression_files/figure-html/unnamed-chunk-14-1.png" width="40%" style="display: block; margin: auto;" /> ] ] ??? * Optionen der `hypothesize()`-Funktion: * "independence" (bivariat) * "point" (univariat, unabhängige Stichprobe (independent measures)) und * "paired independence" (univariat, abhängige Stichproben (repeated measures)) --- class: inverse, center, middle ## Regressionsanalyse im traditionellen Ansatz .blockquote[Koeffizienten und t-Test] .blockquote[Modellgüte und F-Test] .blockquote[Modellauswahl] ??? auf Basis theoretischer Verteilungen --- class: left .blockquote[Koeffizienten und t-Test] ## Regressionstabelle und Ergebnis-Vergleich mit Resampling-Methodik .panelset[ .panel[.panel-name[Anpassung] ``` r library(moderndive) # für get_regression_()`-Funktionen - Ausgabe im tidy-Format # Fit regression model: reg_auto <- tbl_autos %>% lm(Preis~Alter,.) reg_auto %>% get_regression_table() %>% gt() ```
term
estimate
std_error
statistic
p_value
lower_ci
upper_ci
intercept
23056.714
468.871
49.175
0
22107.53
24005.894
Alter
-2635.669
166.935
-15.789
0
-2973.61
-2297.727
] .panel[.panel-name[Vergleich] * Quantitative Übereinstimmung des Koeffizienten `\(b_1=-2635.67\)` mit `obs_slope` aus der **infer**-Pipe * Was ist mit den weiteren Spalten? * Standardfehler der Punktschätzer * t-Wert der Punktschätzer * p-Wert zum t-Test `\((H_0: \beta_0=0\)`, `\(H_A: \beta_1\neq0)\)` * Konfidenzintervalle der Koeffizienten ] .panel[.panel-name[Konfidenzintervalle] <img src="data:image/png;base64,#03_InferenceRegression_files/figure-html/unnamed-chunk-16-1.png" width="55%" style="display: block; margin: auto;" /> ] ] ??? * Anpassung?
**an die Daten** * Hinweis moderndive: * `get_regression_table()`: Tabelle mit Schätzwerten und Statistiken zu den Koeffizienten wie abgebildet * `get_regression_points()`: Tabelle mit Beobachtungen, Vorhergesagten Werten und Residuen * `get_regression_summaries()`: Tabelle mit Statistiken zum Modell (u.a. `\(R^2\)`, F-Wert etc.) * Standardfehler * allg: Standardabweichung der Stichprobenverteilung eines Punktschätzers, dessen Wert für eine Stichprobe berechnet wird. * hier: bezogen auf Koeffizienten --- class: left .blockquote[Koeffizienten und t-Test] ## Standardfehler der Punktschätzer * Allgemein: Standardfehler als Standardabweichung der Stichprobenverteilung eines Punktschätzers * hier also: im Beispiel können wir eine Variation der Steigungsparamter von 166.94 Euro erwarten * Somit: Mass für die Variation der angepassten Koeffizienten
entscheidender Wert für die **Reliabilität** (vgl. t-Test) * Koeffizient `\(b_1=-2635.67\)` (Steigungsparameter) minimiert für die vorliegende Stichprobe die Fehlerquadratsumme. * Für alternative Stichprobe (aus derselben Grundgesamtheit) würde höchstwahrscheinlich ein anderer Wert für `\(b_1\)` die Fehlerquadratsumme minimieren (
**Stichprobenvariabilität**). * Zur Erinnerung: **Bootstrapping**-Verteilung als Approximation dieser Stichprobenverteilung ??? * Standarfehler: im Beispiel können wir eine Variation der Steigungsparamter von 166.94 Euro erwarten Recall from Subsection 8.7.1 that the bootstrap distribution is an approximation to the sampling distribution in that they have a similar shape. Since they have a similar shape, they have similar standard errors. However, **unlike the sampling distribution, the bootstrap distribution is constructed from a single sample**, which is a practice **more aligned with what’s done in real life**. --- class: left .blockquote[Koeffizienten und t-Test] ## Standardfehler und t-Test .panelset[ .panel[.panel-name[Einordnung] *
verwendet den Standardfehler * ...zur Berechnung der Konfidenzintervalle und... * ...für die Durchführung des t-Tests, ob der jeweilige **Koeffizient** sich signifikant von Null unterscheidet: `$$H_0: ~\beta_1=0~~~\mbox{vs.}~~~H_A: ~\beta_1\neq 0$$` * Obs! Zweiseitiger Hypothesentest! * traditioneller vs. randomisierter t-Test ] .panel[.panel-name[Teststatistik] * theoriebasierter (!) Standardfehler des angepassten Koeffizienten `\(b_1\)`: `$$\mbox{SE}_{b_1}=\frac{\frac{s_y}{s_x}\cdot\sqrt{1-r^2}}{\sqrt{n-2}}$$` * Teststatistik `$$t=\frac{b_1-\beta_1}{\mbox{SE}_{b_1}}$$` bzw. unter der Annahme, dass `\(H_0: \beta_1=0\)` wahr ist: `\(t=\frac{b_1}{\mbox{SE}_{b_1}}\)`
ermittelt den p-Wert aus einer t-Verteilung mit `\(n-2\)` Freiheitsgraden. ] .panel[.panel-name[Warum t-verteilt?] * Berechnung Regressionskoeffizient `\(b_1\)` mittels `\(b_1 = \frac{s_{X, Y}}{s^2_X}\)` * Die Kovarianz `\((s_{X,Y})\)` als Linearkombination zweier als normalverteilt angenommener Zufallsvariablen ist dabei normalverteilt (vgl. Annahmen) * Die Varianz `\((s^2_{X})\)` ist unter Annahme normalverteilter Residuen (vgl. Annahmen) `\(\chi^2\)`-verteilt (Summe von quadrierten normalverteilten Zufallsvariablen). * Ein Quotient aus einer standardnormalverteilten Zufallsgröße im Zähler und der Wurzel einer `\(\chi^2\)`-verteilten Größe geteilt durch die Freiheitsgrade im Nenner folgt einer t-Verteilung: `$$t=\frac{b_1-\beta_1}{\mbox{SE}_{b_1}} \sim t(n-2)$$` ] .panel[.panel-name[p-Wert] ``` r reg_auto %>% get_regression_table() %>% select(term,p_value) %>% gt() ```
term
p_value
intercept
0
Alter
0
* p-Wert als Wahrscheinlichkeit für das Auftreten des beobachteten Koeffizienten (-2635.67), wenn die Nullhypothese wahr ist.
im Beispiel ist der p-Wert quasi Null, der Steigungskoeffizient weicht damit mindestens auf dem 1%-Niveau signifikant von Null ab. ] ] ???
Test zur Überprüfung, ob sich der Koeffizient signifikant von Null unterscheidet. Here, our null hypothesis `\(H_0\)` assumes that the population slope `\(\beta_1\)` is 0. If the population slope `\(\beta_1\)` is truly 0, then this is saying that there is no true relationship between teaching and “beauty” scores for all the instructors in our population. In other words, `\(x\)`= "beauty" score would have no associated effect on `\(y\)`=teaching score. The alternative hypothesis `\(H_A\)`, on the other hand, assumes that the population slope `\(\beta_1\)` is not 0, meaning it could be either positive or negative. This suggests either a positive or negative relationship between teaching and "beauty" scores. * t im Auto-Beispiel: -2635.67/166.94 `\(\approx\)` -15.7881275 * Don’t worry if you’re feeling a little overwhelmed at this point. There is a lot of background theory to understand before you can fully make sense of the equations for theory-based methods. That being said, **theory-based methods and simulation-based methods** for constructing confidence intervals and conducting hypothesis tests **often yield consistent results**. vgl. moderndive [Abschnitt 10.5.1](https://moderndive.com/10-inference-for-regression.html#inference-conclusion) * the standard error of $b_1 depends on * `\(\frac{s_x}{s_y}\)`: the relationship between the variability of the response variable and the variability of the explanatory variable . * `\(\sqrt{1-r^2}\)`: Next, it looks into how the two variables relate to each other * Nenner `\(\sqrt{n-2}\)`: most important
as the sample size `\(n\)` increases, the standard error `\(\mbox{SE}_{b_1}\)` decreases * Just **as we demonstrated** in Subsection 7.3.3 when **using shovels** with `\(n = 25\)`, `\(50\)`, and `\(100\)` slots, the **amount of sampling variation of the fitted slope `\(b_1\)` will depend on the sample size `\(n\)`
**Zunehmende Reliabilität (kleinerer Standardfehler) mit zunehmender Stichprobengröße** * First `\(s_x\)` and `\(s_y\)` are the sample standard deviations of the explanatory variable `\(Alter\)` and the response variable `\(Preis\)` respectively. * Second, `\(r\)` is the sample correlation coefficient between `\(Preis\)` and `\(Alter\)`. This was computed as -0.9315177 * Lastly, `\(n\)` is the number of pairs of points in the `tbl_autos` data frame, here: 40. --- class: left .blockquote[Koeffizienten und t-Test] ## Traditionelle Berechnung und Interpretation des Konfidenzintervalls .panelset[ .panel[.panel-name[Vorgehen] * Berechnung mit der Formel: `$$b_1\pm t\cdot \mbox{SE}_{b_1}$$` * Ermitlung des **t-Werts** unter Berücksichtigung... * ...des (grundsätzlich vor der Schätzung der Parameterwerte) Konfidenzniveaus (bspw. 95%) ``` r sig <- .975 #obs! zweiseitiger Test ``` * ...und der Freiheitsgrade (in der Einfachregression: `\(n-2\)`) ``` r df <- pull(count(tbl_autos))-2 # Freiheitsgrade df ``` ``` ## [1] 38 ``` ] .panel[.panel-name[t-Wert in R] ``` r t <- qt(sig,df) t ``` ``` ## [1] 2.024394 ``` ] .panel[.panel-name[Konfidenzintervall in R] ``` r reg_auto %>% get_regression_table() %>% select(term,estimate,std_error) %>% filter(term=="Alter") %>% summarise(ci_lo=estimate-t*std_error, ci_up=estimate+t*std_error) %>% gt() ```
ci_lo
ci_up
-2973.611
-2297.727
Obs! Entspricht (erwartungsgemäss) nicht dem Konfidenzintervall aus dem Bootstrapping-Ansatz (ebenfalls Standardfehler-Ansatz) ] .panel[.panel-name[Interpretation] * **Exakt**: Bei wiederholter Stichprobenziehung (aus der Grundgesamtheit!) können wir erwarten, dass 95% der jeweiligen Konfidenzintervalle den wahren Koeffizienten `\(\beta_1\)` beinhalten. * Unser Konfidenzintervall für den Schätzer `\(b_1\)` (-2973.61; -2297.727)... * ...gehört entweder zu den 95% der Konfidenzintervalle, die den wahren Wert von `\(\beta_1\)` umfassen, sodass dieser nicht Null sein kann. * ...oder es gehört zu den 5% der Konfidenzintervalle, die den wahren Wert nicht umfassen. Dann könnte `\(\beta_1\)` mit geringer Wahrscheinlichkeit auch Null sein. * Kurzform (**nicht exakt**): "Wir sind 95% "zuversichtlich" (*confident*), dass der wahre Koeffizient `\(\beta_1\)` zwischen -2973.61 und -2297.727 liegt." (vgl. ([Ismay and Kim, 2024](https://moderndive.com/index.html)), Abschnitt 10.2.4) ] ] ??? In einem einfachen linearen Regressionsmodell schätzen wir: `\(\beta_0\)` (Achsenabschnitt / Intercept) `\(\beta_1\)` (Steigung des Regressors X) Das sind zwei Parameter, die aus den Daten geschätzt werden. Für jeden geschätzten Parameter geht ein Freiheitsgrad verloren, da wir beim Schätzen die Daten "gebraucht" haben. --- class: left .blockquote[Modellgüte und F-Test] ## Statistiken zur Beurteilung des Modells (Modellgüte) ``` r reg_auto %>% get_regression_summaries() %>% gt() ```
r_squared
adj_r_squared
mse
rmse
sigma
statistic
p_value
df
nobs
0.868
0.864
1884209
1372.665
1408.324
249.281
0
1
40
* (Angepasstes) Bestimmtheitsmass: Misst den Anteil der erklärten Streuung `\((R^2)\)`. `\(\mbox{adj.}~R^2\)` berücksichtigt die Zahl der Regressoren.
Bei starker Abweichung von `\((R^2)\)` und `\(\mbox{adj.}~R^2\)`: Verdacht auf Overfitting * MSE, RMSE und Sigma: kleine Werte bedeuten, dass die Modellvorhersagen gut sind * Statistik und p-Wert: beziehen sich auf den F-Test, der prüft, ob das Modell insgesamt signifikant besser ist als ein Modell ohne Prädiktoren. ??? * Modell ohne Prädiktoren wäre hier einfach der Mittelwert als Schätzer, statt des bedingten Mittelwerts. * df: Freiheitsgrad des Regressors * nobs: Anzahl Beobachtungen Modellzusammenfassung r.squared (Bestimmtheitsmaß R²) 📌 Allgemein: Anteil der Varianz der Zielvariable, der durch das Modell erklärt wird. 🔍 Interpretation: Nur 3.5% der Streuung in mpg kann durch den Prädiktor erklärt werden — sehr schwaches Modell. adj.r.squared (korrigiertes R²) 📌 Allgemein: Wie R², aber angepasst auf die Anzahl der Prädiktoren — schützt vor Überanpassung. 🔍 Interpretation: Bei 3.29 % liegt der bereinigte R² — der Unterschied ist klein, da nur ein Prädiktor verwendet wird. sigma (Residuen-Standardabweichung) 📌 Allgemein: Durchschnittliche Abweichung der tatsächlichen Werte von den vorhergesagten. 🔍 Interpretation: Die durchschnittliche Abweichung liegt bei ca. 0.535 Einheiten auf der Skala von mpg. statistic (F-Wert) 📌 Allgemein: Testet, ob das Modell insgesamt signifikant besser ist als ein Modell ohne Prädiktoren. 🔍 Interpretation: Der F-Wert von 16.7 deutet auf ein signifikantes Modell hin. p.value (F-Test p-Wert) 📌 Allgemein: Wahrscheinlichkeit, dass das beobachtete Ergebnis rein zufällig zustande kommt, wenn kein Effekt vorliegt. 🔍 Interpretation: Mit p = 0.0000508 ist das Modell auf dem 0.1 %-Niveau hochsignifikant. df (Freiheitsgrade für die Regression) 📌 Allgemein: Anzahl der erklärenden Variablen im Modell. 🔍 Interpretation: Es wurde 1 Prädiktor verwendet. logLik (Log-Likelihood) 📌 Allgemein: Maß für die Modellgüte – je höher, desto besser (aber nur im Vergleich mit anderen Modellen sinnvoll). 🔍 Interpretation: Log-Likelihood beträgt –366, allein nicht aussagekräftig. AIC (Akaike Information Criterion) 📌 Allgemein: Modellgüte mit Strafterm für Modellkomplexität – niedrigere Werte sind besser. 🔍 Interpretation: AIC von 738 kann mit anderen Modellen verglichen werden (z. B. mit mehreren Prädiktoren). BIC (Bayesian Information Criterion) 📌 Allgemein: Wie AIC, aber stärkerer Strafterm für Komplexität. 🔍 Interpretation: BIC von 751 ist ebenfalls nur im Modellvergleich nützlich. deviance (Residuenabweichung) 📌 Allgemein: Summe der quadrierten Abweichungen – je kleiner, desto besser passt das Modell. 🔍 Interpretation: Die Residuenabweichung beträgt 132. df.residual (Freiheitsgrade der Residuen) 📌 Allgemein: Anzahl der Beobachtungen minus Anzahl der geschätzten Parameter. 🔍 Interpretation: 461 Freiheitsgrade – deutet auf insgesamt 463 Beobachtungen hin. nobs (Anzahl der Beobachtungen) 📌 Allgemein: Anzahl der in das Modell eingeflossenen Datenpunkte. 🔍 Interpretation: Das Modell basiert auf 463 Beobachtungen. --- class: left .blockquote[Modellgüte und F-Test] ## Vorhersagequalität: MSE und Co. .panelset[ .panel[.panel-name[Vorhersagen & Residuen] ``` r reg_auto %>% get_regression_points() %>% head() %>% gt() ```
ID
Preis
Alter
Preis_hat
residual
1
10000
3.78
13093.886
-3093.886
2
21850
1.61
18813.288
3036.712
3
14500
2.28
17047.389
-2547.389
4
11100
5.33
9008.600
2091.400
5
6700
5.49
8586.893
-1886.893
6
24000
0.34
22160.587
1839.413
] .panel[.panel-name[MSE] * Mittlerer quadratischer Fehler (Mean Squared Error, MSE) als durchschnittlicher quadrierter Vorhersagefehler: `$$\mbox{MSE}=\frac{1}{n}\sum_{i=1}^n(y_i-\hat{y_i})^2$$`
Ein niedriger MSE bedeutet, dass das Modell die tatsächlichen/gemessenen Werte sehr genau vorhersagt. * Problem: MSE in quadratischen Einheiten (hier `\(\mbox{Geldeinheiten}^2\)`) * Lösung: RSME als Wurzel des MSE ] .panel[.panel-name[MSE in R] ``` r mse <- reg_auto %>% get_regression_points() %>% summarise(RSS=sum(residual^2), MSE=1/n()*RSS) mse ``` ``` ## # A tibble: 1 × 2 ## RSS MSE ## <dbl> <dbl> ## 1 75368344. 1884209. ``` ``` r rsme <- sqrt(pull(mse)) rsme ``` ``` ## [1] 1372.665 ``` ] .panel[.panel-name[Sigma] * `\(\sigma\)`: Standardabweichung der Residuen: `$$\sigma=\sqrt{\frac{\mbox{RSS}}{n-p}}$$` ``` r sigma <- sqrt(mse$RSS/(df)) sigma ``` ``` ## [1] 1408.324 ``` ] ] ??? Obs! Werte stimmen mit den Werten aus der `get_regression_summaries()`-Tabelle überein --- class: left .blockquote[Modellgüte und F-Test] ## F-Test auf gemeinsame Signifikanz der Regressoren .panelset[ .panel[.panel-name[F-Test] * Ziel: Überprüfung, ob das Regressionsmodell insgesamt eine signifikante Erklärungskraft hat. * Referenz: Nullmodell (
einfacher arithmetischer Mittelwert) * Hypothesen des F-Tests: - Nullhypothese: Alle Regressionskoeffizienten (außer dem Interzept) sind gleich Null. `$$H_0: \beta_0 = \beta_1 = \dots = \beta_p = 0$$` - Alternativhypothese: Mindestens einer der Regressoren hat einen signifikanten Einfluss. `$$H_A: \text{Mindestens ein } \beta_i \neq 0$$` ] .panel[.panel-name[F-Wert] * Der **F-Wert** wird aus der **Modellgüte** und dem **Fehlermodell** berechnet: - **Modellgüte**: durch das Modell erklärte Variation (Regression). - **Fehlermodell**: nicht durch das Modell erklärte Variation (Fehler). `$$F = \frac{\text{Erklärte Variation (Modell)}}{\text{Nicht erklärte Variation (Fehler)}} = \frac{\frac{SS_{\text{reg}}}{p}}{\frac{SS_{\text{res}}}{n - p - 1}}$$` * Interpretation: - **Höherer F-Wert**: Hinweis auf ein besseres Modell. - **Kleiner F-Wert**: Das Modell erklärt wenig der Variation der Zielgröße. ] .panel[.panel-name[p-Wert] * Vergleich des p-Werts zum F-Test mit dem kritischen Wert: * p-Wert < kritischer Wert bedeutet, dass mindestens ein Regressor einen Einfluss hat... * und das Modell somit einen signifikant höheren Erklärungswert hat als das Nullmodell * Auto-Beispiel: Der p-Wert von 0 deutet darauf hin, dass das Modell mindestens auf dem 1%-Niveau einen signifikanten Erklärungswert gegenüber dem Nullmodell hat. ] ] ??? * Refresher: Regression dient der Ermittlung bedingter Mittelwerte * Referenz des F-Tests also: unbedingter Mittelwert (der abhängigen Variable) Mit anderen Worten: Liege ich beim Schätzen im Durchschnitt besser, wenn ich einfach den Mittelwert (hier: 1.6542\times 10^{4} Euro) der abhängigen Variablen tippe? --- class: left .blockquote[Modellauswahl] ## Schrittweise Selektion .panelset[ .panel[.panel-name[Alternativen] * Rückwärts-Eliminierung: * Beginn mit einem Modell, das alle potenziellen Regressoren enthält * Schrittweise Eliminierung einzelner Regressoren,... * ...bis der Erklärungswert des Modells nicht mehr verbessert werden kann * Vorwärts-Auswahl: * Beginn mit einem Modell, das nur einen Regressor enthält * Schrittweises Hinzufügen einzelner Regressoren,... * ...bis der Erklärungswert des Modells nicht mehr verbessert werden kann * Häufig genutztes Kriterien für den Erklärungswert (vgl. <a name=cite-crh_intro_2024></a>([Çetinkaya-Rundel and Hardin, 2024](https://openintro-ims.netlify.app/))): `\(\mbox{adj.}~R^2\)` als Kriterium * Ausserdem: Kriterium der Sparsamkeit (Parsimony) ] .panel[.panel-name[Beispiel] * Alter und Fahrleistung (km) als Regressoren `$$\begin{eqnarray} \widehat{\mbox{Preis}}&=&b_0 + b_1\cdot \mbox{Alter} + b_2\cdot \mbox{km}\\ \end{eqnarray}$$` ``` r reg_auto2 <- tbl_autos %>% lm(Preis~Alter+km,.) ``` ] .panel[.panel-name[Koeffizienten] ``` r reg_auto2 %>% get_regression_table() %>% gt() ```
term
estimate
std_error
statistic
p_value
lower_ci
upper_ci
intercept
22649.884
411.870
54.993
0
21815.356
23484.412
Alter
-1896.264
235.215
-8.062
0
-2372.854
-1419.674
km
-0.031
0.008
-3.943
0
-0.047
-0.015
] .panel[.panel-name[Modell] ``` r reg_auto %>% get_regression_summaries() %>% gt() ```
r_squared
adj_r_squared
mse
rmse
sigma
statistic
p_value
df
nobs
0.868
0.864
1884209
1372.665
1408.324
249.281
0
1
40
``` r reg_auto2 %>% get_regression_summaries() %>% gt() ```
r_squared
adj_r_squared
mse
rmse
sigma
statistic
p_value
df
nobs
0.907
0.902
1326806
1151.871
1197.658
180.117
0
2
40
Gibt es einen Mehrwert, des erweiterten Modells? Welches Modell sollte gewählt werden? ] ] ??? * Kriterium der Sparsamkeit: Parsimony --- class: left .blockquote[Modellauswahl] ## Multikolinearität * In Modellen mit mehreren Regressoren kann das Problem der Multikollinearität auftreten. * Multikollinearität als (nahezu) lineare Abhängigkeit zwischen zwei oder mehr Regressoren. * Folge: größere Standardfehler, unsichere Schätzungen (aber keine systematische Verzerrung!) | Problem | Führt zu verzerrten Schätzern? | Ursache | |----------------------------------|-------------------------------|-------------------------------------------| | **Endogenität / Identifikationsproblem** | ✅ Ja | Korrelation zwischen Regressor und Fehlerterm | | **Multikollinearität** | ❌ Nein (aber hohe Varianz) | Lineare Abhängigkeit zwischen Regressoren | ??? Multikollinearität betrifft nicht die Gültigkeit der Inferenz im engeren Sinne, sondern ihre Zuverlässigkeit Die genannten vier Annahmen (Linearität, Unabhängigkeit, Normalverteilung, Homoskedastizität) sichern ab, dass die Schätzer unverzerrt und die Tests korrekt sind. Multikollinearität verletzt keine dieser vier Annahmen direkt. * Aber: Sie führt dazu, dass Standardfehler der Koeffizienten groß werden, was zu unsicheren Schätzungen und breiten Konfidenzintervallen führt. * Deshalb ist Inferenz formal noch gültig, aber praktisch oft wenig aussagekräftig. --- class: left .blockquote[Modellauswahl] ## Überprüfung auf Multikolinearität .panelset[ .panel[.panel-name[Korrelation] ``` r tbl_autos %>% summarise(r_XY=cor(Alter,km)) ``` ``` ## # A tibble: 1 × 1 ## r_XY ## <dbl> ## 1 0.797 ```
sehr hoher Wert (Anfangsverdacht, aber keine hinreichende Bedingung!) ] .panel[.panel-name[VIF] * Varianzinflation (Variance Inflation Factor, VIF) ``` r car::vif(reg_auto2) ``` ``` ## Alter km ## 2.745211 2.745211 ``` * Faustregel: * VIF < 5: unproblematisch * VIF > 5: Vorsicht geboten * VIF > 10: Problematisch ] ] ??? * Obwohl Alter und Kilometerstand korreliert sind (cor = 0.80), bleibt die Multikollinearität beherrschbar. Die Schätzer bleiben stabil und interpretierbar. * Bei problematischen VIF-Werten: bspw. Lösung durch **Hauptkomponentenanalyse** und Verwendung der Hauptkomponenten als Regressoren (
Verweis auf ER017) --- class: inverse, center, middle ## Regressionsdiagnostik im traditionellen Ansatz .blockquote[Annahmen für valide Inferenz] .blockquote[Ausblick: BLUE-Eigenschaft] --- class: left .blockquote[Annahmen für valide Inferenz] ## Regressionsdiagnostik zur Überprüfung kritischer Annahmen .panelset[ .panel[.panel-name[Annahmen] * Linearität: Es besteht ein linearer Zusammenhang zwischen Prädiktoren und der Zielvariablen. * Unabhängigkeit: Residuen sind nicht systematisch voneinander abhängig (z.B. bei Zeitreihen). * Normalverteilung: Die Residuen sind (annähernd) normalverteilt. * Homoskedastizität: Die Varianz der Residuen ist über alle Werte von `\(\hat{y}\)` konstant. ] .panel[.panel-name[Visualisierung] * Mit dem **ggfortify**-Paket können vier Plots zur Regressionsdiagnostik in einer Zeile erzeugt werden: <img src="data:image/png;base64,#03_InferenceRegression_files/figure-html/unnamed-chunk-33-1.png" width="35%" style="display: block; margin: auto;" /> * Interpretationen zu den Plots auf den folgenden Folien ] ] --- class: left .blockquote[Annahmen für valide Inferenz] ## Überprüfung der Linearitäts-Annahme .panelset[ .panel[.panel-name[Beispiel] ``` r tbl_autos %>% ggplot(aes(x=Alter,y=Preis)) + geom_point(alpha=.5,size=.8) + theme_light() + geom_smooth(method = "lm", se = FALSE) ``` <img src="data:image/png;base64,#03_InferenceRegression_files/figure-html/unnamed-chunk-34-1.png" width="35%" style="display: block; margin: auto;" /> ] .panel[.panel-name[Visuelle Inspektion] ``` r library(ggfortify) autoplot(reg_auto2, which = 1) ``` <img src="data:image/png;base64,#03_InferenceRegression_files/figure-html/unnamed-chunk-35-1.png" width="30%" style="display: block; margin: auto;" /> * Interpretation: * Punkte sollten zufällig um die Nulllinie streuen
Linearität erfüllt * Systematische Muster
Hinweis auf nichtlineare Zusammenhänge ] ] ??? Interpretation der visuellen Inspektion: * Streudiagramm: Residuen vs. Vorhersagewerte * Ziel: Kein systematisches Muster, zufällige Streuung um die Nulllinie * Kurven oder Muster deuten auf Verletzung der Linearität hin --- class: left .blockquote[Annahmen für valide Inferenz] ## Beispiel für Verletzung der Linearitäts-Annahme <img src="data:image/png;base64,#https://moderndive.com/ModernDive_files/figure-html/non-linear-1.png" width="75%" style="display: block; margin: auto;" /> .quelle[Quelle: ([Ismay and Kim, 2024](https://moderndive.com/index.html))] --- class: left .blockquote[Annahmen für valide Inferenz] ## Unabhängigkeit der Residuen .panelset[ .panel[.panel-name[Visuell] * Autokorrelationsplot der Residuen <img src="data:image/png;base64,#03_InferenceRegression_files/figure-html/unnamed-chunk-37-1.png" width="35%" style="display: block; margin: auto;" /> * Interpretation: systematische Muster in den Residuen gibt, die auf eine Verletzung der Unabhängigkeit hinweisen ] .panel[.panel-name[Test] * Durbin-Watson-Test für Autokorrelation ``` r library(lmtest) dwtest(reg_auto) ``` ``` ## ## Durbin-Watson test ## ## data: reg_auto ## DW = 2.8695, p-value = 0.9982 ## alternative hypothesis: true autocorrelation is greater than 0 ``` Interpretation: Ein Wert nahe 2 deutet darauf hin, dass keine Autokorrelation vorliegt. Werte deutlich unter oder über 2 deuten auf Abhängigkeiten hin. ] ] ??? * Durbin-Watson: `\(H_0\)` (keine Autokorrelation)... * ...kann beim hohen p-Wert = 0.9982 nicht abgelehnt werden kann.
Dies bedeutet, dass keine signifikante Autokorrelation der Residuen existiert. --- class: left .blockquote[Annahmen für valide Inferenz] ## Normalverteilung der Residuen .panelset[ .panel[.panel-name[Visuell] * Quantile-Quantile-Plot (QQ-Plot) <img src="data:image/png;base64,#03_InferenceRegression_files/figure-html/unnamed-chunk-39-1.png" width="35%" style="display: block; margin: auto;" /> Interpretation: Abweichungen vom Diagonalstrich deuten auf Abweichungen von der Normalverteilung hin. ] .panel[.panel-name[Test] * Shapiro-Wilk-Test auf Normalverteilung ``` r shapiro.test(residuals(reg_auto)) ``` ``` ## ## Shapiro-Wilk normality test ## ## data: residuals(reg_auto) ## W = 0.97287, p-value = 0.4416 ``` * Shapiro-Wilk `\(H_0\)` (Normalverteilung der Residuen)... * kann beim hohen p-Wert = 0.44 nicht abgelehnt werden kann. ] ] ??? Interpretation: Der Q-Q-Plot zeigt, ob die Residuen ungefähr einer Normalverteilung folgen. Abweichungen vom Diagonalstrich deuten auf Abweichungen von der Normalverteilung hin. W = 0.97287: Keine signifikante Abweichung von der Normalverteilung. p-Wert = 0.4416: Nullhypothese (Normalverteilung) wird nicht abgelehnt. Ergebnis: Residuen sind normalverteilt und erfüllen die Annahme der Normalverteilung. --- class: left .blockquote[Annahmen für valide Inferenz] ## Homoskedastizität .panelset[ .panel[.panel-name[Visuell] * Residuen vs. Vorhersage <img src="data:image/png;base64,#03_InferenceRegression_files/figure-html/unnamed-chunk-41-1.png" width="35%" style="display: block; margin: auto;" /> Interpretation: Wenn die Residuen keine erkennbare Struktur (bspw. Trichterform) aufweisen , spricht dies für Homoskedastizität. ] .panel[.panel-name[Test] * Breusch-Pagan-Test für Homoskedastizität ``` r bptest(reg_auto) ``` ``` ## ## studentized Breusch-Pagan test ## ## data: reg_auto ## BP = 1.1565, df = 1, p-value = 0.2822 ``` * Breusch-Pagan `\(H_0\)` (Homoskedastizität liegt vor)... * kann beim hohen p-Wert = 0.28 nicht abgelehnt werden kann. ] ] ??? * BP = 1.1565: Teststatistik. * p-Wert = 0.2822: Keine Heteroskedastizität festgestellt. * Ergebnis: Die Annahme der Homoskedastizität ist erfüllt. --- class: left ## Exkurs: Ausreisser-Diagnostik * Der vierte Plot aus `autoplot(reg_auto)` (sog. Leverage-Plot) wurde bisher nicht genutzt. <img src="data:image/png;base64,#03_InferenceRegression_files/figure-html/unnamed-chunk-43-1.png" width="35%" style="display: block; margin: auto;" /> * Beobachtungen, die außerhalb der gestrichelten Linien (die Schwellenwerte für Einfluss) liegen, haben einen großen Einfluss auf das Modell und sollten überprüft werden. * Hohe Werte in Cook's Distance (Ordinate) oder Leverage zeigen potenziell problematische Punkte, die die Modellgüte verfälschen könnten. --- class: inverse,center,middle ## Ausblick .blockquote[Multiple Regression und BLUE-Eigenschaften] .blockquote[Logistische Regression] --- class: left .blockquote[Multiple Regression und BLUE-Eigenschaften] ## BLUE-Eigenschaft (Best Linear Unbiased Estimator) * In der multiplen Regression kann Multikollinearität auftreten und zu instabilen Schätzungen und erhöhten Standardfehlern führen. * Unter den Annahmen der klassischen linearen Regression (Linearität, Unabhängigkeit der Fehler, Homoskedastizität, Normalverteilung der Fehler) sind OLS-Schätzer BLUE (**B**est, **L**inear **U**nbiased **E**stimator)
OLS-Schätzer sind unter den o.a. Annahmen bestmögliche unverzerrte Schätzer `\((E[\hat{\beta}] = \beta )\)` mit den geringsten Standardfehlern. * vgl. dazu ausführlich <a name=cite-vonAuer_oekonometrie_2023></a>([von Auer, 2023](https://link.springer.com/book/10.1007/978-3-658-42700-9)) ??? (und damit höchste Präzision innerhalb der Klasse der validen Schätzer) --- class: left .blockquote[Logistische Regression] ## Binärer statt numerischer Regressand * In der Praxis interessiert häufig die Vorhersage binärer Ergebnisse: * bspw. Klassifikation von E-Mails in "Spam" oder "Nicht-Spam" * Vorhersage von Kaufentscheidungen ("Kauf" oder "Nicht-Kauf") * Bei binären abhängigen Variablen (z.B. Ja/Nein) kann die lineare Regression nicht angewendet werden, da die Vorhersage außerhalb des Intervalls [0, 1] liegen kann. * Lösung: logistische Regression auf Basis einer logistischen Funktion --- class: inverse,center,middle # Schönen Feierabend. --- background-image: url("data:image/png;base64,#http://bit.ly/cs631-donkey") background-size: 80% --- class: left ## Quellenverzeichnis .ref-slide[ <a name=bib-vonAuer_oekonometrie_2023></a>[Auer, L. von](#cite-vonAuer_oekonometrie_2023) (2023). _Ökonometrie. Eine Einführung_. Wiesbaden: SpringerGabler. ISBN: 978-3-658-42699-6. URL: [https://link.springer.com/book/10.1007/978-3-658-42700-9](https://link.springer.com/book/10.1007/978-3-658-42700-9). <a name=bib-crh_intro_2024></a>[Çetinkaya-Rundel, M. and J. Hardin](#cite-crh_intro_2024) (2024). _Introduction to Modern Statistics_. OpenIntro. ISBN: 978-1943450275. URL: [https://openintro-ims.netlify.app/](https://openintro-ims.netlify.app/). <a name=bib-ismay_statistical_2024></a>[Ismay, C. and A. Y. Kim](#cite-ismay_statistical_2024) (2024). _Statistical Inference via Data Science: A ModernDive into R and the Tidyverse_. 2nd ed. The R Series. Boca Raton London New York: CRC Press, Taylor & Francis Group. URL: [https://moderndive.com/index.html](https://moderndive.com/index.html). ]